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A two-component system of penetrable particles interacting via a gaussian core potential is con- 
, sidered, which may serve as a crude model for binary polymer solutions. The pair structure and 

. thermodynamic properties are calculated within the random phase approximation (RPA) and the 

' hypernetted chain (HNC) integral equation. The analytical RPA predictions are in semi-quantitative 

agreement with the numerical solutions of the HNC approximation, which itself is very accurate for 
gaussian core systems. A fluid-fluid phase separation is predicted to occur for a broad range of 
potential parameters. The pair structure exhibits a nontrivial clustering behaviour of the minority 
component. Similiar conclusions hold for the related model of parabolic core mixtures, which is 
frequently used in dissipative particle dynamics (DPD) simulations. 
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I. INTRODUCTION 



Demixing of binary or multicomponcnt mixtures is a very common phenomenon observed in a broad range of 
molecular fluids polymer solutions and blends [§[ ||, or colloidal dispersions Q ||. Phase separation is generally 
associated with differences in the attractive interactions between particles of different chemical species. In polymer 
^ \ solutions these differences are usually embodied in the Flory ^-parameter , which controls the competition between 
Q 1 the entropy of mixing and the total interaction energy, at least within a mean-field picture. On the other hand, in multi- 
O component colloidal systems like binary dispersions involving colloidal particles of very different sizes, or mixtures of 
colloidal particles and non-adsorbing polymer, phase separation can be driven by purely repulsive, excluded volume 
I ' interactions. By mapping the initial multi-component system onto an effective one-component system involving only 
^ , the bigger colloidal particles, the largely entropy-driven demixing can be understood in terms of attractive depletion 
' interactions induced between the large particles by the smaller species (the "depletant" ) |^ . It should however be kept 
00 in mind that the initial bare interactions are purely repulsive, albeit strongly non-additive, as in the highly simplified 
Asakura-Oosawa model for colloid-polymer mixtures In the case of fully additive hard sphere mixtures, phase 
separation has been predicted for sufficiently large size ratios but it is now generally believed that the fluid-fluid 
demixing is metastable, and preempted by freezing A significant degree of positive non-additivity of the core radii 
R^i, (whereby R12 = {Rn + i?22)(l + A)/2, with A > 0) is required to observe a stable demixing transition in the 
fiuid phase 



d Effective interactions between the centres of mass of fractal objects, like linear polymer coils lid, 11, [13, [13], star 



polymers pM or dendrimers |15 , obtained by averaging over individual monomer degrees of freedom, are now known 



I to be very "soft" . More specifically the effective pair potential diverges only logarithmically for overlapping star 



■ O polymers [ |l^ , while remaining finite, of the order of 1 — 2kBT, for linear polymers in good solvent | |lO| , 



IIL 



This observation has stimulated the investigation of simple models, like finite repulsive step potentials or the 
gaussian core potential ||l^, |l^, |l9j, which was first introduced by Stillinger, in a somewhat different context [ pO[ , 
namely 

w(r) = eexp(-rVi?^) , (1) 



o 
o 

> 
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■ where e is the energy scale, while R determines the range of the effective potential. It is worth stressing that the 
"gaussian core" model is unrelated to the "gaussian molecule" model, which was extensively studied by Michael Fisher 
and collaborators In the latter model it is the Mayer /-function, rather than the pair potential, which has a 

gaussian shape. 

Simple, penetrable particle models are also widely used in highly coarse-grained simulations of large-scale phenomena 
within the so-called "dissipative particle dynamics" (DPD) method [E2l|2^. In DPD, effective interactions between 
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penetrable fluid "particles" are frequently modelled by a simple parabolic potential |2J| : 




It has recently been realized that binary mixtures of particles with penetrable cores, which interact via generalizations 
of the gaussian and parabolic potentials and (^), involving different energy scales e and radii R for the various 
species, may phase-separate over appropriate ranges of these parameters. Spinodal instability was first shown to 
occur for the gaussian core model (|l|) within the random phase approximation (RPA) and binodals as well as 
intcrfacial properties were then calculated within the same approximation . Similarly Gibbs ensemble Monte Carlo 
simulations have very recently shown that binary systems of soft particles interacting via the parabolic potential (|^) 
phase separate beyond a critical degree of enhanced repulsion between particles of different species, in agreement with 
Flory-like mean field considerations |Q . 

In this paper we systematically extend our earlier results for the gaussian core model [l^ and investigate the 
range of validity of the RPA by detailed calculations of the pair structure, thermodynamics, and the resulting phase 
coexistence curve within the much more accurate hypernetted chain (HNC) approximation. The break-down of the 
RPA is quantified in the physically relevant regime where pR^ ~ 1 and e ~ fc^T, which would correspond to the 
cross-over from dilute to semi-dilute regimes of the underlying binary polymer solution. The RPA continues to provide 
reliable first estimates at higher densities. 



II. RPA AND HNC 

The model under consideration is the binary gaussian core model (GCM) already introduced in references fl^ and 
[ p^ . It consists of A^i particles of "radius" Ri and N2 particles of "radius" R2 in a volume V. The total number 
density is p = (iVi -l-iV2) /V, while the concentrations of the two species are x = N2/N and 1 — x — Ni/N, respectively. 
The pair potentials are 

^^M.W = ep.e-('-/«-)', (3) 

which introduce three length and three energy parameters: Rn, R22, R12, en, £12 and £22- If ^11 is chosen as unit of 
length, the system is entirely specified by the 5 dimensionless parameters R12/R11, R22IR11, e*! — f^sii, 6*2 — P^i2i 
£22 — P^22, where /3 = \/{kBT). For a fixed set of dimensionless parameters the reduced Helmholtz free energy per 
particle, / — F / {NkBT), is a function only of the intensive variables p and x; this may be split into the ideal gas , 
ideal mixing and excess (non-ideal) contributions: 

f{p, x) /id(p) + /mix (a;) + /ox(/3, x) (4a) 
= In(pA^) ~l + x\xix+{\-x) ln(l ~ x) + fc^{x,p), (4b) 

where A is an irrelevant de Broglie thermal wavelength. The equation of state f3P/ p, and the chemical potentials 
/iy, = 1,2) follow from the standard thermodynamic relations: 



P \dp 



Pf^^ = (¥) --(¥-) (5b) 



\9pJ.^ \dxj 



\9pJ.^ \dx 



p 



and may likewise be split into ideal and excess contributions. Within RPA and HNC, these excess contributions may 
be easily expressed in terms of the usual total and direct correlation functions h^i, (r) and c^^ (r) , which are related by 
the familiar Ornstein-Zernike (OZ) relations j2^; these may be used to express the Fourier transforms (FT) h^y{k) 
of the total correlation functions in terms of the FT c^i,{k) of the direct correlation functions: 

hii{k) = ^ [cii(fc)(l - P2C22{k)) + P2C?2(fc)] , (6a) 
'hi2{k) = ^^ci2(fc), (6b) 
h22{k) = ^ [c22(fc)(l - Picii(fc)) + Pic?2(fc)] , (6c) 
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where pi = p(l — x), p2 = px, and 

A(fc) = [1 - picn(fc)][l - p2C22(fc)] - PiP2cL(fc)- (7) 

These OZ relations must be supplemented by closure relations. 

The RPA amounts to simply identifying the Cf^^{r) with their asymptotic (large r) behaviour: 

= -^^.e-C^/^-)^ (8) 

Substitution of the FT's of (||) into the relations (||) and (^ yields explicit expressions for the hu^{k), which may be 
transformed back to obtain the (r) . Detailed analytical expressions are given in Appendix |A| in the special case 
of mixtures with equal core radii Rn ~ R22 = R12 = R- 
The HNC closure relations are |^ : 

9f,,^{r) = 1 + h^„{r) (9a) 

= exp{-(3v^^{r) +7^^(r)} 
Ifj^i^ir) = hf^^{r) - Cf_,^{r) (9b) 

An iterative procedure must be used to solve the coupled equations (|^) and (^) numerically. Convergence is easily 
achieved for the binary GCM model over most of the potential parameter space. In fact the HNC approximation 
becomes exact in the high density limit (pi?^ — > 00), and is extremely accurate for densities pR^ « 1, as shown earlier 
in the one-component case [|l^ |l8|. Pair structure data will be examined in greater detail in section 5. 

Knowledge of the pair correlation functions allows direct access to thermodynamics, via the compressibility or virial 
routes Being approximate, the RPA and HNC closures are not thermodynamically consistent, i. e. the two routes 
lead to different answers. In the case of the GCM the two theories become strictly thermodynamically consistent only 
in the high density hmit, ^ 

Within the RPA, the compressibility route is equivalent to a simple mean- field ansatz for the free ener gy IT^ , p^ . 
This leads immediately to the following analytic expression for the excess part of the reduced free energy [|lq: 

fe.{^,p)^lpVo{^) (10a) 



lpJ2'^^i'^''^t^'' ^^^^^ 



V fJ, 

V^, = f3v^Ak = 0) = / P^t^Ar)dr, (10c) 
where xi = (1 — x), X2 = x. In the special case of the GCM: 

V^. = (11) 

The pressure and chemical potentials p^ follow directly from eqns. (^ (the superscript C refers to the compressibility 
route) : 

P ^ 

(3p^ = ln(pA^(l - x)) + Viip(l -x) + Vi2px (12b) 
The virial route is considerably more arduous. The pressure follows from the virial theorem: 



= InipA^x) + Vi2p{l -x) + V22PX (12c) 



l + ^EE^A^^-V (13) 
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In the case of the symmetric mixture {Rn = R22 = R12 = R', Vn = V22) the RPA pair correlation functions obtained 
in Appendix A may be used to yield the following expression (the superscript V referring to the virial route): 



/3P^ _ (3P^ Z2x{l - x)A 



p 7r3/2i?3 
zix{l — x)A 



P_ 



^3/2^3 



Z2 



(14) 



where the functions H(s) and the roots z\ and Z2 are defined in appendix^ (eqns. (A4a) and ( A4h| ) ) . The resulting 
reduced free energy per particle is obtained by integrating the equation of state (14) with respect to density: 



r(x,p)=/C(x,p) + 



[Li5/2(p/2:i) + Li5/2(p/^2) + pVi] 



(15) 



where the polylogarithm Li5/2 is defined in equation (A7). Unfortunately this expression becomes singular for densities 
larger than the spinodal densities obtained by the compressibility route. This mathematical artifact discussed in 
Appendix ^ does not allow us to use expressions (|l^) and (^5|) to construct a phase diagram. 

The lack of thermodynamic consistency of the RPA is illustrated in Figure 0, which compares the compressibility 
and virial equations of state as a function of density, for a symmetric mixture; the virial pressures are always lower 
than their compressibility counterparts, and closer to the nearly exact HNC results, to which we now turn. 

The HNC approximation is very nearly thermodynamically consistent in practice, at least for temperatures and 
densities relevant for polymer solutions [ |l8| . Hence we have calculated thermodynamic properties within the more 
convenient virial route with the pressure given by the standard relation (p^, while the chemical potentials may also 
be directly expressed in terms of the direct and total correlation functions, according to |2q]: 



- fir(p,A3) + I '^'^^"^ 



(16) 



Note that these expressions hold only within the HNC approximation, and are consistent with the virial route (O 



III. SCALING PROPERTIES 



As mentioned earlier, suitably reduced equilibrium properties of the binary GCM depend on the five dimcnsionless 
combinations R12/R11, R22/ Rii, e*!, 6*27 ^^^d €22- The compressibility version of the RPA thermodynamics (or equiv- 
alently, the mean field approximation) allows a considerable reduction of this parameter space. From the expressions 
(p^)-(12c) it is clear that the radii P^j, and energies e^i^ enter only in the combinations V^j^y defined in eqn. (11). In 
terms of the reduced density p* = pVn and pressure P* = pPVn, the thermodynamic behaviour of the mixture is 
uniquely determined by the dimcnsionless ratios V22/V^ii and Vi2/^/Vn\^. This remarkable reduction of parameter 
space holds only for the thermodynamics, but not for the correlation functions, which depend explicitly on the five 
parameters of the GCM potentials. The scale invariance of the thermodynamics is broken by the RPA virial route, as 
is immediately evident form the explicit expressions (14) and ( [T^ ) (valid for a symmetric mixture); the same is true 
within the HNC approximation. 

However deviations from RPA-C (mean-field) scale invariance are expected to be small in the high density, high 
temperature regime, where RPA is increasingly accurate, so that the considerable reduction in the number of relevant 
potential parameters (from 5 to 2) is expected to carry over, at least approximately, to the phase diagrams calculated 
with RPA-V or HNC thermodynamics, which will be presented in the next section. 



IV. PHASE DIAGRAMS 



The phase behaviour of the binary GCM may be deduced from the knowledge of the reduced free energy per 
particle / as a function of the variables v = 1/p and x and from its thermodynamic derivatives (the pressure and the 
chemical potentials). For suitable values of the potential parameters e^^^ and the binary GCM becomes unstable 
against demixing at sufficiently high densities. In terms of the free energy per particle, / — f(v,x), the standard 
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thermodynamic stability conditions of a binary mixture read Q : 



^'^^ >0; (^) >0 (17a) 

V 

2 

> 0. (17b) 



\ '^'"^ J x\ J V \ dvdx J 



The first condition ensures mechanical stability, while the second inequality guarantees stability against demixing at 
constant volume; the third inequality ensures stability of the mixture at constant pressure. Note that the stability at 
constant pressure is the more restrictive condition. 

The subsequent discussion is restricted to demixing at constant pressure. The vanishing of the 1. h. s. of equation 



(171), corresponding to the case where the third inequality turns into an equality, signals the occurrence of spinodal 
instability. Within the mean- field (RPA-C) approximation (|l^), the equation for the spinodal is easily calculated to 
be 

1 + pVi {x) ~ p^x{l- x)AV === 0, (18) 

where: 

Fi(a;) = (1-x) 1^11 +2:^22 (19a) 

Ay = y22-yiiy22. (i9b) 

It is easily inferred from (^8|) that a spinodal instability occurs whenever Vi2/V^ii'^22 > 1- This region is visualized 
in Figure || in a plot of the minimum energy ratio ei2 / Y^eiie22 required for spinodal instability, as a function of the 
size ratio i?22/-Rii- The radius R12 is taken to be given by the combination rule: 

^?2 = ^(^?i+^22), (20) 

as suggested by renormalization group (RG) calculations ||ll|] and by direct simulation ||lj] of mixtures of self-avoiding 
polymer coils. 

Figure || also shows the values of the ratio ei2/-^eiie22, as calculated for two self-avoiding walk polymers by RG 
techniques |l^, plotted as a function of the ratio of their radii of gyration. The figure suggests that phase separation 
of real polymers would only be observed for ratios of gyration radii larger than 10, i. e. for mixtures of very long and 
very short polymers. However, the RG results are probably not trustworthy for such asymmetric mixtures, and the 
RG potentials are strictly valid only in the infinite dilution limit. 

The spinodal line in the P — x-plane is easily calculated from (p^) to be: 



Vi{x) + ^[V,{xW+4x{l-x)AV 
^ 2x{l-x)AV 



The critical concentration is determined by the condition: 

dPsjx) 
dx 



= 0, (22) 



where Ps{x) is the pressure calculated from (12a), at the spinodal density given by (pT|). The corresponding critical 
density pc, calculated by substituting the critical concentration determined by ( |2^ ) is plotted in Figure ^ versus the 
ratio Vi2/V^ii^22, for several ratios V22/V11. 

The binodal or phase coexistence curves in the p ~ x plane are determined by the usual conditi ons o f equa lity of 



the chemical potentials of both species and of the pressures in the two phases, using expressions ( 12a )-(12c). The 
resulting equations must in general be solved numerically to yield the binodal curves. The complete phase diagram can 
however be calculated analytically within the RPA-C approximation, in the symmetric case, where Vn — V22. Due 
to the symmetry of the problem, all thermodynamic quantities must be invariant with respect to the transformation 
1 <-* 2; X ^ {1 — x) and p ^ p. This simplification allows the fully analytic treatment detailed in Appendix B. The 
resulting spinodal and binodal are shown in Figure]^ in the case where yi2/Vii ~ 2. The same Figure also shows the 
HNC binodals, for several values of eli (remember that within RPA-C all these binodals coincide). As expected the 
HNC results break the scale invariance of RPA-C and shift the binodal curves upward, to higher densities. The HNC 
and RPA-C critical densities differ by almost a factor of 2 for eli = 2, pointing to the limitations of the mean-field 
(RPA-C) description. 
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An example of a phase diagram in an asymmetric mixture, with potential parameters already used in earlier RPA 
calculations jl^, |l9|, is shown in Figure || in the p — x plane and in Figure || in the P — x plane. The spinodal and 
binodal curves are now asymmetric, and the difference between RPA-C and HNC results is again significant, with the 
HNC binodal being pushed to higher densities. 

Phase separation in binary mixtures of particles interacting via the closely related parabolic core potentials has 
been observed in the Gibbs ensemble Monte Carlo simulations of Wijmans et al. j24|. They considered the symmetric 
case where all i?^^ are equal, = and e*2 = e^^i + Ae*. We have calculated the binodal in the Ae* —x plane under 
the same conditions as the MC simulations, i. e. for = 12.5, = pbR^ = 3, in the RPA and HNC approximations. 
The results are compared to the MC data in Figure |7[ In view of the fact that the amplitude of the parabolic repulsion 
is 12.5^5^1 which makes it more hard-core like, the agreement may be considered to be rather satisfactory. As in the 
GCM case, the HNC coexistence-curve lies well above its RPA counterpart; it is closer to the simulation data. 



V. PAIR STRUCTURE AND CLUSTERING 

The pair correlation functions hfj,,y{r) of the GCM mixture can be calculated by combining the OZ relations (^) with 
either the RPA closure (||) or the HNC closure (^). In the symmetric case, characterized by a single gaussian range 
parameter Rn = R22 — R12 — R, analytic expressions, in the form of infinite series, can be obtained for the RPA 
closure, as shown in Appendix A. In the asymmetric case, solutions for the ft.pi/(r) are readily obtained by numerical 
Fourier transformation of equations (^) and , using the straightforward Fourier transforms of the c^^ (r) given by 

The HNC closure requires an iterative solution of the coupled OZ and closure equations. This was achieved using 
the standard Picard method, and well-converged solutions were generally obtained in a few iterations. Convergence 
was found to be slower in the vicinity of phase coexistence, and to break down rapidly inside the phase coexistence 
region corresponding to metastable mixtures. Concentration fluctuations build up in that region, as signalled by a 
k = peak of growing amplitude in the h^^(k), or equivalently in the corresponding partial structure factors 

phf,.{k) (23) 

The concentration-concentration structure factor: 

Scc{k) = xlS^m - 2xiX2Si2{k) + xlS22ik) (24) 
satisfies the long wavelength limit |25[ : 

lim S'cc(fc) = tr.2r^it^I\ ' (25) 

fc^o {d^G/dx^)N.p,T 

where G is the Gibbs free energy. The cross-over from metastable to unstable mixture corresponds to the vanishing 
of {d'^G/ dx^)N,p,T, so that Scc{k = 0) is expected to diverge along a spinodal hue. 

This is indeed the case with the RPA closure, but the HNC closure ceases to converge before a spinodal line is 
reached, a well-known shortcoming of the HNC closure [ p7[ . Such a break-down is typical of the inadequacy of fluid 
integral equations to describe critical behaviour ||2^. The accuracy of HNC in the stable one-phase region for the 
pair correlation function /i(r), already documented in the one-component GCM ||l^, |l^, is tested under more severe 
conditions (namely e* = 12.5) in Figure || against recent Monte Carlo data |Q for the parabolic core potential (H). 
The agreement is seen to be excellent. This suggests that the HNC closure would be very useful to generate pair 
structure in coarse-grained DPD fluids at a moderate computational cost. 

Examples of pair correlation functions and partial structure factors for several states of GCM mixtures are shown 
in Figures p[]Tl|. The most striking feature is the amount of structure observed both in r and in /c-space, even within 
the RPA, compared to the previously studied one-component case iQ. The peak in the structure factors at fc = as 
the spinodal is approached is expected, as explained earlier. More surprising perhaps is the appearance of peaks in the 
pair correlation functions as r — ^ 0, or at finite r. Peaks at r ^ are of course precluded in the presence of hard cores, 
but are quite significant in the present penetrable core model and point to a novel physical mechanism to trigger phase 
separation, namely the clustering of the minority species, as observed e. g. in Figure ^for the symmetric case. The 
clustering can be explained by a simple "energetic" argument: the particles of the majority species tend to maximise 
their mutual distances. If the repulsion between unlike particles is stronger than between like ones, particles of the 
minority species prefer to cluster in the voids left by the majority species, rather than overlap with particles of the 
latter. This mechanism drives phase separation. 
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VI. CONCLUSION 

The two-component extension of the gaussian core model leads to non-trivial phase behaviour, which we have 
investigated within RPA and HNC theories. The former is partially analytic, while the latter requires only a modest 
numerical effort to calculate partial pair correlation functions and the resulting thermodynamic properties. While the 
pair structure shows some unexpected features, the most interesting prediction is the occurrence of phase separation 
induced by purely repulsive pair interactions. As expected for this model of penetrable particles, phase separation is 
driven by an enhanced repulsion between unlike particles compared to the like-particle interaction. 

A special feature of the GCM mixtures is that at high densities of both species, the mixture behaves like a "mean 
field fluid", i. e. the RPA becomes asymptotically exact jl^, [l8|. This means that at finite densities, the RPA makes 
semi-quantitavely valid predictions for the phase behaviour, and allows a rapid exploration of potential parameter 
space to search for likely conditions for phase separation. 

The nature of the link between the behaviour of the binary GCM and that of a binary polymer solution must 
still be worked out in detail, mainly because the effective interactions between the centres of mass of the polymer 
coils are state-dependent |L3| . For a one-component system at lower densities, the structure is not very sensitive to 
the state-dependence [p9|, pO|. We therefore expect the binary GCM to make qualitatively correct predictions for 
the structure of binary polymer solutions in the dilute and beginning of the semi-dilute regimes. The link to the 
phase-behaviour is less clear. Because the GCM uses a fixed potential, it is mainly relevant to behaviour in the dilute 
regime. While polymers in a melt are known to phase-separate rather easily |Q, in the dilute regime for good solvents 
they are not expected to phase-separate [||. This appears to be confirmed by the results shown in Fig 2. However, it 
may be that for poorer solvent, or for stronger polymer incompatibilities, a regime where the polymers phase-separate 
at low densities could open up. If this is the case, then the binary GCM would be a useful coarse-grained model with 
which one could rapidly explore the qualitative phase behaviour of this regime. 

Another open question relates to the critical behaviour of the binary GCM in the immediate vicinity of the critical 
consolute point. The RPA clearly predicts mean field exponents. As regards the exact exponents, it is not obvious 
whether the correct critical behaviour of the model necessarily belongs to the Ising universality class, since there is no 
clear correspondence between GCM mixtures and standard lattice models, due to the penetrability of the particles. 
The question of the correct universality class, a subject dear to the heart of Michael Fisher, must be considered as 
open. 
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APPENDIX A: SOLUTION OF THE VIRIAL ROUTE IN RPA 



We consider the GCM where the potentials have equal range Rn = R12 = R22 = R, and the heights of the repulsion 
potentials of like particles are the same for both species en — £22- Without loss of generality we can set Vu = 1. 
Then we have 

c^.(fc) = -V;.-e~'■'^'/^ (Al) 

with 

V^.. = (A2) 

If we introduce the abbreviation z{k) = pV\ie~^ ^ 1^ and set AV^ = Vj^j ~ ^ii — ^12 ^ 1j the OZ equations (^) reduce 
to (using equation (^): 

f X2AVz(kY ~ zik) , 

Phii{k) - ^ (A3a) 

P^i2(fc) = -^ (A3b) 

^ x^AVzjkf - z{k) 

ph22(k) = -^-^ (A3c) 

A{k) = 1 + z{k) - x{l- x)AVz{kf (A3d) 

The denominator A(fc) in equation (||) has the zeros 



^1 = Vt; Vttt^ > (A4a) 

2x{l-x)AV ^ ' 

Z2 = — -, < (A4b) 

2x{l--x)AV ^ ' 



The positive root z\ is the spinodal density at any given composition. 
The virial equation ( p^ gives 



27r /■ 
= — -^P^^^i^^^ \ '"'^"l^^'^' ' h^^{r)dr (A5a) 

= J^EE^M^"/ rfre^.(r-V^')e"'''^'''VW- (A5b) 



Evaluating this integral in Fourier space and making use of the Parseval theorem leads to 

"3 A:2i?2 

/J- 1/ 



2 4 



= — / [a:ia;2£(fc) AV^ - 2xia;2 A^ - 1] [6 - k^R^] (A6) 

12 J A(fc) 

This integral cannot be evaluated analytically in general. However, if the density is below the spinodal density for 
a given composition, one can expand the denominator in the integrand in a power series with respect to z{k). The 
integral can then be evaluated term by term, interchanging the integration and the summation. In order to make the 
symmetry between the two species explicit, we write a; — The shifted concentration ^ now varies between — 1 

and 1, with ^ = corresponding to x — 0.5. 

With the definition of the polylogarithmic function 
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and the auxiliary function 



Ha) = — [Li3/2(-a) - Li5/2(-a)] 



evaluation of ( [A^ ) leads to 



P 



1 + -p[i + Vv2 - i'Av] 



47r3/2i?3 

4^3/2^3 

This can be integrated to give the free energy per particle 

Pf{x,p) = f'' + r-{x) 



p_ 

Z2 



1 



[Li5/2(p/2l) + Li5/2(/?/z2) + p] 



From that we get the chemical potentials as 

ppi = log(A3p) + log 



+ 



1 



+ 



27r3/2i?3(i _ 

1 



T- ( p 
L13/2 ( - 



27r3/2i?3(i_^)^' 
fip2 = log(A3p) + log 



1+n , 1 + e 



+ 



Li3/2 I — 



27r3/2i?3(i + ^)^ 

2^3/2^3(1 + L'3/2 



+ V^12 


1 + e 


1 


1 


2 ^ 


27r3/2i?3 












- Ay(i 








- Ay(i 






+ V^12 




1 


1 


2 ^ 


27r3/2i?3 












- Ay(i 


-a. 






- Ay(i 







(A8) 



(A9) 



(AlO) 



(All) 



(A12) 



It is easily checked that = M2(— due to the symmetry of the model. Note that these expressions 

become singul ar w hen p> z\ = p^. This is a mathematical artefact originating in the interchange of summation and 
integration in (A6), which renders the expressions for P and p unusable for binodal calculations. 
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APPENDIX B: ANALYTIC SOLUTION OF THE RPA-C BINODAL IN THE SYMMETRIC CASE 



The phase diagram simpUfies considerably in the symmetric case where Vu — ¥22- We express V12 as V12 — V\\ (I+7) 
and note that demixing occurs when the dimensionless parameter 7 is positive. As above we introduce the reduced 
quantities p* = pV\\ and P* = (iPV\\. Because of the symmetry of the problem, all quantities must be invariant 
with respect to the transformation p — > p, a; — > 1 — x, 1 <-> 2, which corresponds to a relabeling of species 1 and 2. 
As in appendix |a|, we make use of the shifted concentration 1^ to make this symmetry explicit. The pressure and the 
chemical potentials can now be expressed in terms of the variables p* and ^ : 

P*=p* + ip*2(2 + 7)- (Bla) 
= ln(p*) + In i^^f) + ^P*(2 + 7) + \p*il (Bib) 
/3p2=ln(p*)+ln(^i±^^ +ip*(2 + 7)-^P*^7 (Blc) 

In these equations we have set A — Vii, which merely shifts the free energy per particle by a constant amount and 
thus does not change the phase diagram. 

Two coexisting phases a, /? must have the same pressure and chemical potentials, respectively: 

= P\ = p^ = p^2 (B2) 

Because of the symmetry in this model, coexisting phases must have the same density p" = p^ and the corresponding 
concentrations must fulfill x" — 1 — . Considering this, equal pressure in both phases is trivially fulfilled, and 
the two conditions of equal chemical potentials become equivalent. The relation between the binodal density pi, and 
concentration is 

P% = 4ln(i±|). (B3) 
Expanding the binodal densities around the critical point gives 

(B4) 



2 

7 



1 + + Oi^) 



From the series expansions we see that for all 7 > the spinodal lies inside the binodal. Within the RPA approximation 
we thus get an analytically exact solution for a fluid-fluid demixing transition of a liquid consisting of purely repulsive 
particles. As expected in this treatment, the resulting critical exponents are mean-field exponents. 
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FIG. 1: The normalised effective second viri al coefficient B2 = 4(/3P - p)/(3pVii) is plotted versus the normalized density 
p*/2 at fixed concentration x — 0.5 for RPA-C, RPA-V and HNC . In the upper plot, we have e*i = 622 = 0.1, e*2 = 0.2, while 
in the lower plot eji = £22 = 1-0, £12 = 2.0. Both systems have the reduced critical density p* = 2 in RPA-C. The RPA-C 
route predicts B2 = 1 at all densities (solid line). The RPA-V values are represented by the dotted lines, the HNC data by 
the dashed lines up to the density beyond which convergence fails. The agreement between the three routes to the equation of 
state is seen to improve with increasing density, except in the vincinity of the critical density where the results diverge again. 
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FIG. 2: Mixing/demixing regions within RPA-C, in the ei2 / y/^\\^22 vs. R22/R11 plane. Above the solid line phase separation 
is predicted by RPA-C. Corresponding values of ei^ A/ei 1 622 predicted for polymers with radii of gyration _Rii and R22 by 
a renormalisation group theory of Kriiger et al. |ll| are plotted as the dashed curve. The intersection of these curves at 
R22/R11 ~ 10 suggests the possibility of a phase separation in polymer mixtures with extreme size ratios. 
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FIG. 3: Reduced critical densities PcVu vs. Vii/VVuT^, for several values of ^22/^11 from RPA-C (shown in the inset). 
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FIG. 4: Binodal for a symmetric GCM mixture with Vi2/\^22 = 2 for values of en = £22 ~ 0.2,0.5,1.0,2.0, determined by 
RPA-C (full curves) and HNC (dashed curve). The RPA-C binodal and spinodal (dotted line) depend only on the ratio Fi2/Vii. 
The HNC binodal does not obey the same scale invariance and is shifted to higher reduced density as en increases. 
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FIG. 5: Phase diagram in p — x plane for asymmetric mixture with = 1, R\2 = 0.665, R22 = 0.849, £11 = £22 = 2, £*2 = 1.888. 
The sohd curve is the RPA-C binodal with examples of tie-lines of coexisting phases as dash-dotted lines. The dotted line is 
the RPA spinodal. The circles and dashed curve represent the HNC binodal with tie-lines shown again as dash-dotted lines. 
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FIG. 7: Binodal of the parabolic core system (g) investigated by Wijmans et al. [gj]. Here en = £22 ~ 12.5 is kept fixed 
and e\2 is determined so that the binodal density has the fixed value pbR'' = 3.0 at any given concentration x. The required 
difference Ae* — ~ ^ii is plotted vs. x. Note that in reference twice that value is plotted due to a different definition of 
the potential parameters. The diamonds correspond to the simulation data of Wijmansei al. [^^, the solid line represents the 
RPA-C prediction. The HNC results are plotted as circles. While RPA-C underestimates the required for the given binodal 
density by a factor of two, the HNC results lie closer to the simulation data, although a significant discrepancy remains. This 
may be due to the high amplitude £^2 of the potentials, making the particles more hard-sphere like. 
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FIG. 8: A comparison between HNC predictions for gi2{r) and the simulation data of Wijmans et. al. for the parabohc 
potential with eJi = = 12.5. The potential amplitude between unlike partikles e\2 is varied as 612 = en + Ae* with 
Ae* = 0, 6.25, 12.5. The simulation data are plotted as dashed lines, while HNC results are plotted as solid lines. Note that 
(712 ^ 0) is very near zero, in contrast to the case of the much softer potentials studied in llq] and the subsequent Figures. 
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FIG. 10: HNC partial structure factors for the nori-syrnriietric system Rn = 1.0, i?i2 ~ 1.582, i?22 = 2.0, e*i = ei2 = £22 = 
2.0, pRii = 0.188, X = 0.5. The solid lines represents the RPA partial structure factors, the dashed lines the HNC results. The 
enhanced values at fc = show that this system is close to the spinodal. 



20 




r/R 

FIG. 11: HNC pair correlation functions of particles in the symmetric case with e^i = £22 = 0.1, £^2 = 0.2, at a concentration 
X = 0.9 and densities below the binodal {pR^ = 2.0, solid line), on the binodal {pR^ = 5.0, dotted line) and near the spinodal 
{pR^ = 6.7, dashed line). The results for the minority species 1 point towards clustering. 



